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OVER BLUNT-NOSED AND FLARED BODIES 

By Mamoru Irxou ye, John V. Rakich, and Harvard Lomax 
Ames Research Center 


SUMMARY 


The computer programs developed at Ames Research Center for calculating 
the inviscid flow field around blunt-nosed bodies are described briefly and 
their application to specific shapes is demonstrated. The programs solve 
numerically the exact equations of motion for plane or axisymmetric bodies at 
zero angle of attack and for a perfect gas or a real gas in thermodynamic 
equilibrium. An inverse method is used for the subsonic-transonic region, and 
the method of characteristics is used for the supersonic region. Results are 
shown for several body shapes in both perfect and real gas flow, including a 
comparison between air and a C0 2 -N 2 mixture . Presented are shock-wave shapes 
and distributions of pressure and other flow variables along the body and 
across the shock layer. 


INTRODUCTION 


Aircraft and spacecraft designers are faced with the problem of determin- 
ing the inviscid flow field over blunt-nosed bodies for supersonic flight at 
speeds encompassing those attained in planetary entry. In addition to the 
blunt nose, a typical body shape may have a flared afterbody which further 
complicates the problem. The dominant features of such a flow field are indi- 
cated in sketch (a) . There occurs a detached bow wave that is normal at the 



Sketch (a) 



axis of symmetry and decays in strength as it approaches a Mach wave at large 
distances from the body. The flow behind the shock wave is subsonic in the 
nose region bounded by the sonic line and becomes supersonic over the after- 
body. Expansion waves and embedded shocks may occur as a result of corners. 

In addition, an embedded shock may arise from coalescence of compression waves 
from the surface or from separation of the boundary layer, which often occurs 
on this type of body. Analysis of the viscous region is beyond the scope of 
the present study; however, its study depends on a knowledge of the external 
inviscid flow. 

A number of exact and approximate techniques for determining the flow 
field depicted in sketch (a) have been reported. Some of the more recent con- 
tributions are references 1 through k- for blunt-body flows, and references 5 
and 6 for the supersonic region downstream of the nose. For flared bodies, 
exact numerical results have been reported in reference 7 while approximate 
methods may be found in references 8 through 10. Hayes and Probstein (ref. ll) 
present a more complete discussion of the entire flow field. 

The computer programs that are described in the present report solve 
numerically the exact equations of motion for plane and axisymmetric flow at 
zero angle of attack and provide the complete inviscid flow field between the 
body and the shock wave . The fluid may be a perfect gas or a real gas in 
thermodynamic equilibrium. An inverse method (ref. 3) is used for the 
subsonic-transonic region (referred to as the blunt-body program) , and the 
method of characteristics (ref. 5) is used to extend the calculations down- 
stream in the supersonic region. These computer programs were written in 
FORTRAN II for use on an IBM 709^ a t Ames Research Center, but have been made 
available to a number of other organizations. The distribution of these pro- 
grams has created a need for a more complete description and documentation 
than is presently available. The present report is intended to partially 
fulfill this need. 

The purpose of the present report is to provide a general description of 
the Ames flow-field computer programs and to present results of calculations 
that demonstrate the range of applicability. The governing equations of 
motion are introduced briefly at the start. Then the methods used to solve 
the equations are presented. No attempt is made to provide a complete listing 
of all the subroutines and flow charts. Instead, detailed descriptions are 
provided only for selected portions of the programs that warrant special con- 
sideration. The information contained in this report should acquaint the 
reader with the general logic followed in the programs and be helpful in diag- 
nosing small difficulties or in making minor modifications. 

Sample results are presented for shock-wave shape, surface- pres sure 
distribution, and shock-layer profiles of total pressure, static pressure, 
density, and velocity for various free-stream conditions and body shapes. 

The first examples demonstrate how a simple modification improves the accuracy 
of the calculations in regions with large entropy (or vorticity) gradients. 

Then comparisons are made with flow-field results obtained by means of an 
integral method for the blunt-body solution. Comparisons are also made with 
experimental results obtained for a body with a flare. Finally, examples of 
calculations for real gases in thermodynamic equilibrium are presented. 
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SYMBOLS 


a speed of sound 

ellipsoid bluntness, (b/c) 2 
b,c semiaxes of ellipsoid 

h enthalpy 

total enthalpy 
M Mach number 

n coordinate normal to a streamline 

p static pressure 

p^ total pressure 

R nose or cylinder radius 

S entropy 

s,t sheared coordinates (see eqs . (8)) 
u,v velocity components in x,y directions 
V velocity 

X sho c k- wave shap e 

x,y cylindrical coordinates with origin at body nose 

7 ratio of specific heats 

A shock standoff distance 

5 angle of corner on the body 

€ index for number of degrees of symmetry; e = 0 for plane symmetric 

flow, and e = 1 for axisymmetric flow 

0 flow angle 

0 C cone angle 

M- 


Mach angle 



p density 

i|r stream function 


Subscripts 


b body 
s shock 

oo free-stream conditions 


EQUATIONS 


The partial differential equations that must be satisfied for steady, 
inviscid flow are as follows: 

Continuity of mass 

■1 (puyt + -1 (pvy e ) = 0 
ox by 

where e = 0 for plane symmetric flow and e - 1 for axisymmetric flow. 
Momentum equations 
x direction 


(i) 


+ P v — + 
dx by dx 


= 0 


y direction 


dv , „„ dv , dp n 
pu — + pv — + — — = 0 


dx 


dy dy 


Energy equation 


( 2 ) 


(3) 


bx by 


U + v _ a 2( u ^P + v SR) = 0 


dp d 

— + v - 
dx dy 


where a is the isentropic speed of sound defined by 

v<W s 


a 2 = fa 




(5) 


k 



To solve these equations for a given set of boundary conditions, the 
thermodynamic properties of the gas are required; for example, 

a 2 = f(p,p) (6) 

For a perfect gas these relationships are merely functions of the gas constant 
and ratio of specific heats* For example, equation (6) becomes 

*2 = f ( 7 ) 


For a real gas the equilibrium composition and thermodynamic properties must 
be obtained by means of statistical mechanics. These calculations can be done 
independently of the flow-field equations, and the results can be tabulated 
for later use. Dr. Harry E. Bailey of Ames Research Center has recently per- 
formed these calculations for various gas mixtures of current interest follow- 
ing the assumptions and approximations made by Marrone (ref. 12). The data 
cover temperatures to 25,000° K in 250° increments and densities from 10“ T to 
10 s times a reference density, p Q , which is the density of the mixture for a 
temperature of 273*l6° K and a pressure of 0.101325 MN/m 2 (l atmosphere). 

For example, the properties for carbon dioxide are reported in reference 13* 

The thermodynamic properties in this form are not suitable for optimal 
use in a computer program. Some approximations are necessary to minimize the 
computing time and storage requirements. For use in the present programs, the 
calculated values of the properties have been spline fitted with cubics by the 
method of reference l4, and the coefficients of the cubics have been stored on 
magnetic tape. A special subroutine reads the tape, searches for the proper 
coefficients, and evaluates the desired properties. This approximate tech- 
nique, in general, yields results within 1 percent of the original data. At 
present the thermodynamic properties for air and the twelve mixtures of nitro- 
gen, carbon dioxide, and argon listed in table I are available on tape. 

For moderate temperatures, for example, below about 2000° K for air, 
dissociation and ionization can be neglected, and the imperfect gas effects 
are due to the excitation of the vibrational states. The thermodynamic prop- 
erties for such thermally perfect gases have been calculated in reference 15 
and have also been stored on tape for use in the present programs . 

The system of equations is now complete. In general, the four partial 
differential equations (eqs. (l) through ( 4 )) must be solved simultaneously 
for the four dependent variables p,p,u, and v. In the following sections 
the methods used to solve these equations numerically in the subsonic- 
transonic region and in the supersonic region are discussed. 


METHOD OF SOLUTION FOR SUB SONIC -TRANSONIC REGION 


In the nose region of blunt bodies, equations (l) through (^) exhibit 
different character; namely, the equations are elliptic in the subsonic region. 
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parabolic on the sonic line, and hyperbolic in the supersonic region. Despite 
these complications, an inverse method (see, e.g., ref. l) has been found 
effectual for solving such flow fields. In the application of this method, a 
shock shape is assumed and the equations are integrated numerically by a 
finite-difference method to determine the corresponding body shape. The par- 
ticular version used in this report is reported in detail in reference 3; 
hence, only a brief description will follow. 


Since the initial boundary conditions are specified along the shock, a 
sheared, nonorthogonal coordinate system with one axis coincident with the 
shock is useful (see sketch (b)) . The new coordinates are defined as follows: 



s = x - X(y) 

t = y 


(8) 


Equations (l) through (4) are then 
transformed and expressed in the following 
form: 


ap 

5s 

5s 


5u 

5s 




>V, 


p,u,v, 



p,u,v, 



p,u,v. 





p,u,v. 


dp 

ap 

du 

dv\ 

ah 

at' 

at' 

at/ 

dp 

dp 

du 

dv\ 

dt' 

at' 

at' 

at/ 

dp 

dp 

du 

8 v\ 

ah 

at' 

at' 

at/ 

dp 

dp 

du 

a v \ 

a-h 

at' 

at' 

at/ 


(9) 


For a given set of free-stream conditions and shock shape, the values of 
p,p,u, and v/t just behind the shock wave are calculated from the Rankine- 
Hugoniot relations, and the derivatives with respect to t are found by 
numerical differentiation. Then equations ( 9 ) are used to march in step-by- 
step toward the body. A flow chart for the computer program is showi in 
figure 1, and a list of subroutines is provided in table II. To illustrate 
the predictor-corrector integration technique, suppose the flow properties 
are known for the (i - l)th and ith steps and are to be calculated for the 
(i + l)th step (see sketch (b)). A second-order predictor and a modified 
Euler ian second-order corrector are used as follows, where p is a typical 
flow variable . 


1. Differentiate numerically data for ith step to obtain (5p/5t)^ 

2. Calculate from equations ( 9 ), (5p/5s)^ 
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3 . Predict new value p 1 = p. + 2 As(3p/ds). 
b. Differentiate numerically to obtain (dp/dt)^ +i 
5 • Calculate from equations (9), (dp/ds)^ +i 
6. Correct value p. , = p. + 0.5 As[(dp/ds). + (dp/ds) ] 

-Li-L X X X+X 

The stream function is calculated for each pointy and the body is deter- 
mined as the locus of points where the stream function vanishes. The step 
size is chosen so that the stagnation point is reached in approximately seven 
steps. It is usually necessary to iterate on the shock shape to obtain the 
desired body shape. However, this procedure is simplified by a one-parameter 
family of shock shapes that will produce reasonably accurate spheres and 
ellipsoids . Values of the shock-shape parameter for spheres are presented in 
reference 3 for perfect gases and for air in thermodynamic equilibrium. Solu- 
tions for other gases and nonspherical bodies including not-too-blunt ellip- 
soids can be obtained by appropriate changes of the shock-shape parameter. 
There are limitations on the application of the program, mainly because of 
inherent numerical difficulties, as discussed in reference 3. As a general 
rule, solutions are not possible for Mach numbers less than 3 and for body 
shapes that are either very blunt (B 5 > b) or not smooth in the nose region. 

The output from the blunt -body program consists of the flow properties on 
the body and in the flow field at the coordinate intersections shown in 
sketch (b) . In addition, the properties are interpolated along a line joining 
the shock and body in the supersonic region. These- data can be used as input 
for the method of characteristics program to continue the calculations down- 
stream over the afterbody. 


METHOD OF SOLUTION FOR SUPERSONIC REGION 


A computer program based on the method of characteristics is used to 
determine the flow field in the supersonic region. This program is comprised 
of a main program and 33 subroutines which are listed in table III. Most of 
these subroutines are short and straightforward and, therefore, they will not 
be explained in detail. (Table III describes the primary function of each.) 
However, the main program requires a few words of explanation. In addition, 
certain quadratic interpolation procedures as well as methods for calculating 
embedded shock waves will be discussed. 


Calculation Procedure for Smooth Bodies 

Along the characteristic or Mach lines, the partial differential equa- 
tions ((l) through ( 3 )) reduce to the following ordinary equations (see, e.g., 
ref. ll): 

g sin 0 gy /j_q\ 

m sin(e ± p) y 


c°|Ji dp ± ae = - 
pV 2 
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Equation (10) is solved in conjunction with the energy equation in integrated 
form 

+ — = h-fc = constant (ll) 


h 


the conservation of entropy on streamlines 

S = S(n) 


( 12 ) 


and the equation of state 


h = h(p,S) 


( 13 ) 


Shock 


in array, P (J,K ), 
where 

J - identities flow 
variables (x,y, V,— etc) 
K- identifies the location 
of the field points 
(numbered as shown) 

Direction of calculation 


Starting 
line 


Briefly, the method consists of starting with flow properties along a non- 
characteristic line Between the Body and the shock wave, as determined from 
the solution for the suBsonie-transonic region, and then integrating the 
equations downstream along the Mach lines. The stepwise procedure is illus- 
trated By the typical characteristic mesh shown in sketch (c) . Beginning with 

known data on the starting line 
Current data is stored the calculation proceeds to the 

Body along a right-running 
characteristic, and then Back 
to the next starting point (or 
shock point) . A simple flow 
chart is shown in figure 2 
which illustrates this part of 
the program logic . In sketch 
sketch (c), the previously cal- 
culated (or input) data points 
are identified By small circles, 
and the point currently Being 
calculated is identified By the 
shaded symbol. Only the num- 
bered points are available in 
computer memory at this time, 
the remaining circled points 
having Been written out pre- 
viously. The stored data points are contained in a two-dimensional array, 
P(J,K), in which the index. J identifies the various flow variables and the 
index K identifies the location of the point. In terms of program termin- 
ology, the number of field points involved in the calculation loop is given 
By M2 + 1, where M2 is an integer defined in the main program. The field 
point currently Being calculated is identified as P(j,K9)j where K9 is also 
an integer defined in the main program. Thus as the calculation proceeds 
along a right -running characteristic, the integer K9 takes on successive 
values Between K9 = M2 + 1 on the shock (or input line) to K9 = 1 on the 
Body. 

For the calculation of a typical mesh point (C in sketch (c)), three 
adjacent points are usually used. These points are labeled A,D,B in 
sketch (c), and correspond, in the example shown, to the points K = 2,3,^ in 



Typical characteristic mesh 

Sketch (c) 
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the P array. The calculation of data at the new point is effected with the 
use of equations ( 10 ) through (13) and a standard predictor-corrector proce- 
dure which averages the coefficients of the differentials. The procedure is 
started with a crude predictor (i.e., that conditions at C equal those at 
B) and is followed, therefore, by at least two correctors. This is in con- 
trast to the method used for the blunt -body solution which makes use of only 
one corrector, but which uses a second-order predictor. 

In calculating rotational supersonic flow by the method of characteris- 
tics, it is convenient to introduce entropy as a flow property since it 
remains constant on streamlines. In the past it has sometimes been assumed 
that the entropy varies linearly between streamlines (see ref. 16, p. 636). 

To illustrate this procedure, consider three points in the flow field (see 
sketch (c)) A and B, where the flow properties are known, and C, where the 
flow properties are to be determined. The entropy at C can be calculated 
using the flow properties at A and B and with the assumption that the 
entropy varies linearly along the normal to the streamlines between A and B. 
This assumption is valid provided that the change in entropy gradient between 
A and B is small. However, serious errors may occur in the flow-field cal- 
culations if this condition is not realized (see, e.g., ref. 17 )- Decreasing 
the mesh size by increasing the number of mesh points is not a satisfactory 
remedy because the computing time increases as the square of the number of 
input points, and the storage capacity of computing machines is also limited. 
In reference 17, an iterative scheme was used wherein a check on the inte- 
grated mass flow was made point by point throughout the field. The method 
adopted in the present program is simpler in that no additional iterations are 
needed. Quadratic rather than linear interpolation for the entropy is used, 
and errors are therefore reduced to the order of the cube of the mesh size. 
This is accomplished by using the flow properties at point D, which lies on 
Mach lines upstream from A and B. This additional point allows the use of a 
quadratic calculation for entropy between A and B. Some improvement in the 
flow-field solution can be expected, especially in regions where large changes 
in entropy gradient occur. 

A feature of the present computer program, which proved useful for exam- 
ining entropy gradients, is the ability to interpolate for field data on 
radial (body- or axis-normal) lines. This operation is noted in the flow 
chart (fig. 2 ). These interpolated data at several x stations are printed 
out at the end of the normal output and, if desired, the data from the last 
x station can be stored on magnetic tape. This tape can later be used to 
provide starting conditions to extend the calculation downstream. 


Calculation Procedure for Bodies With Corners 

The method of characteristics described in the previous section is not 
applicable in regions where the body slope is discontinuous, or where embedded 
shocks occur. For these cases, characteristics theory may be applied on both 
sides of the discontinuity with matching conditions obtained for the Rankine- 
Hugoniot shock relations, or Prandtl-Meyer equations in the case of an expan- 
sion corner. Methods for calculating flows of this type may be found, for 


iu 
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example, in reference l6. This section describes some of the details of this 
calculation as incorporated into the present program. 

While any number of discontinuities are allowable in theory, practical 
considerations have resulted in a limitation, for the present program, of any 
two of the discontinuities indicated in sketch (a) . These are: (a) an expan- 

sion corner, (b) a compression corner, and (c) an embedded shock arising, for 
example, from a coalescence of Mach waves from a concave wall. Procedures 
used in calculating these three types of discontinuities are discussed next. 
Descriptions of the compression corner and coalesced shock explain only the 
starting conditions used for the shock calculation, and are followed by an 
explanation of a general embedded shock point. The calculation for an embed- 
ded shock proceeds to its intersection with the bow shock. Then, since inter- 
actions between shocks are not considered, the calculation terminates along a 
right-running characteristic through this point. 

Expansion corner .- The expansion corner is illustrated in sketch (d) . 

Upon reaching the body in region II for the first time, point C is calculated 

on the extension of the body shape 
specified for region I. This pro- 
vides an analytic continuation of 
the flow ahead of the corner and, 
with stored data at points A and 
B, enables one to use a quadratic 
interpolation for conditions at 
point D just ahead of the corner. 
Now, given conditions at D and 
the expansion angle 5, the 
Prandtl-Meyer equations can be 
used to calculate conditions on 
the body just behind the corner. 

In addition, conditions for sev- 
eral intermediate angles are com- 
puted and all are stored with 
coordinate values corresponding to 
point D. The problem is now 
reduced to one which can be handled by the main characteristics program. With 
known conditions at points D and E, point F can be computed, followed by 
similar calculations for points G,H, etc., until the entire expansion fan is 
determined as shown in the inset. A greater number of mesh points are intro- 
duced at such a corner for large expansion angles, so as to provide a reason- 
ably uniform mesh. 

Compression corner .- The compression corner is shown in sketch (e) . The 
procedure for calculating conditions at point D is identical with that 
described for the expansion corner. In this case, the oblique shock relations 
are used to calculate the flow variables on the body just behind the corner in 
terms of upstream conditions at point D and the known deflection angle, 5. 

The necessary shock solution is not explicit, however, and an iterative pro- 
cedure has been programmed to give the jump conditions . The segment of the 
shock D-F is assumed straight and at an angle corresponding to the oblique 
shock at point D. Data at point F are then found by a quadratic 
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Shock 



Compression corner 


Sketch (e) 

interpolation along the right-running characteristic through E ; and the jump 
conditions at point F are computed from the shock relations. The data at 
D f and F f can now be used to determine the body point C T . Knowledge of the 
shock wave at point F must be stored in the main computer program so that 
the general shock-point subroutine can be called when the calculation along 
the next right -running characteristic reaches this point. 


Coalesced shock .- In sketch (f) 
flow field is depicted at point E 



General shock point 

Sketch (f) 


the formation of a coalesced shock in the 
where two Mach lines of the same family 
have intersected. A shock wave is 
started at point E and at an angle 
equal to the average slope of the Mach 
lines F-E and G-E. The jump condi- 
tions at A-A* are then computed and 
stored in the P array. Now the points 
numbered (K9 - 3 ) , (K9 - *0 > and so on, 
are computed with the use of known data 
on the characteristic through point G. 
This procedure gives the starting con- 
ditions for the coalesced shock, and 
the problem is now reduced to the gen- 
eral case which is explained next . 


General shock point .- This sub- 
routine solves, in an iterative proce- 
dure, for the shock angle at point B 
in sketch (f ) , which satisfies the flow 


conditions behind the shock. The procedure is as follows. When the shock is 
reached, point C is first calculated in the usual way with data at (K9 - l) , 
K9, and (K9 + l)- After an initial guess, the average shock angle at A and B 
is used to locate point B, and data there are obtained by quadratic inter- 
polation through points C, (K9 + l ) 9 and (K9 + 2). Data at B* are obtained 
from the shock relations, and the intersection D of Mach lines through A 1 
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and B* is determined. Data at point D are obtained by interpolating through 
points A*,(K 9 -3 ), and (K9 *- 4). Equation (10) is then used to determine 
conditions at B' in terms of data at D, that is, on the line D-B' 

£2* d P + de = € sln -g — ^ 

pV 2 M sin(0 + p) y 

The value of pressure at B* obtained from this equation is then compared 
with that obtained from the shock conditions and the iteration is continued 
until the two values agree. The calculation on the downstream side of the 
shock can then proceed in the usual way. 


RESULTS AND DISCUSSION 


Effect of Entropy Calculation on Solutions for a Perfect Gas 

The flow field around a blunted cone is quite sensitive to the entropy 
calculation because of the growth of the body shape. Since the total amount 
of high entropy fluid remains constant, the thickness of the entropy layer 
must diminish as the body radius grows . Consequently, the entropy gradients 
normal to streamlines must increase. 

The flow field over a spherically blunted cone was calculated by the 
present method for 7 = 1.4 and = 10. Both the shock-wave shape and 
surface-pressure distribution were unaffected by the entropy calculation. 
However, the total pressure profiles in the shock layer are significantly 
different as shown in figure 3(a)* (The ordinate is the distance from the 
surface normalized by the total distance to the shock.) The differences are 
associated with the overexpansion of the shock wave which occurs near 
x/R = 10. The entropy rise across the shock wave is a minimum there, and, 
hence, the total pressure in the shock layer has the maximum value of 
2 .lX10 4 Xp o3 on the streamline passing through this point. At every downstream 
station, the total pressure on this streamline must have the same value. The 
linear entropy interpolation smooths the total pressure distribution and 
erroneously reduces this maximum value. Past x/R = 40, the entropy layer 
becomes so thin that even quadratic interpolation is not accurate near the 
body. This deficiency is usually unimportant from a practical standpoint 
because the inviscid entropy layer would have been absorbed within the 
boundary layer. v 

The static pressure profiles are not affected by the entropy calculation 
as shown in figure 3(t>). However, the entropy calculation has a significant 
effect on the density and velocity profiles shown in figures 3(c) and 3(d) . 

The result is that the mass-flow balance between the shock layer and the free 
stream deteriorates with distance as shown in figure 4. At x/R =25, the 
error in the mass-flow balance is 9 percent with linear interpolation and 
negligible with quadratic interpolation. 
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Similar comparisons have been made for a hemisphere cylinder for which 
the entropy gradients are not as large- Linear interpolation for the entropy 
yielded satisfactory results for the shock-wave shape , surface -pres sure dis- 
tribution, shock-layer profiles, and mass-flow balance. These results cor- 
roborate the findings of reference 5 where linear interpolation for the 
entropy was used, and checks of the mass-flow balance showed satisfactory 
agreement. These examples show that for most bodies linear interpolation for 
the entropy is adequate for obtaining shock-wave shape and static -pres sure 
distribution, but it may cause serious errors in the other flow-field proper- 
ties in regions with large entropy gradients. Since quadratic interpolation 
does not require an appreciable increase in the computing time, it is used 
exclusively in the present method. 


Comparison of Results From Present Method 
With Chushkin and Shulishnina 

It is of interest to compare results from the present method with other 
calculations. Flow-field calculations by the method of characteristics for 
blunted cones by Chushkin and Shulishnina (ref. 6 ) were recently brought to 
our attention. These calculations differ from the present method in that the 
subsonic-transonic solution was that obtained by Belotserkovskii (ref. 4) 
using the method of integral relations. In addition, some differences are to 
be expected in the computational procedure for the method of characteristics. 
Solutions were obtained for spherically blunted cones with half angles ranging 
from 0° to 40° and for y = 1.4 and = 4, 6 , and 10. The surface -pres sure 
distributions in the nose region from reference 6 and the present method are 
in good agreement as shown in figure 5« Although not shown in the figure, 
this good agreement extends downstream until the sharp-cone values are reached. 

Comparisons of pressure distributions over cylinders with ellipsoidal 
noses are shown in figure 6 for 7 = 1.4 and k^ = 6 . For the slender or pro- 
late ellipsoid, the calculations proceeded smoothly, and the results from the 
present method show good agreement with the results of reference 6 . For the 
blunt or oblate ellipsoid, the blunt-body solution used in the present method 
is in error in the transonic region. However, the method of characteristics 
solution quickly corrects this error, at least, as far as the surface pressure 
is concerned . 


Comparison With Experiment for a Body With Corners 

As an illustration of the embedded shock and expansion fan calculations, 
the flow over a blunt cone-cylinder-flare body was calculated for perfect air 
at a Mach number of 4.10. The results of the shock shapes obtained from this 
calculation are shown in figure 7, and the surface pressures are shown in 
figure 8 . Also shown in these figures are experimental results from refer- 
ence l 8 . Of particular interest in figure 7 is the good agreement between 
experiment and theory for the embedded shock on the flare . Good agreement is 
also found for the experimental and calculated pressure distributions in 
figure 8 . Note that the calculation predicts an increasing pressure with dis- 
tance along the flare. This pressure variation is caused by the increase of 
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upstream dynamic pressure with distance from the corner, that is, the effect 
of the nonuniform upstream conditions created by the strong bow shock (cf. 
ref. 9) . 


Example of Solution for Gas Mixtures 
in Thermodynamic Equilibrium 

Flow-field solutions for gases in thermodynamic equilibrium present no 
inherent difficulty aside from an increase in computing time required to cal- 
culate the thermodynamic properties. However, the limited accuracy of the 
real gas properties obtained from the curve fits coupled with the different 
calculation methods used in the blunt body and method of characteristics pro- 
grams may result in an incompatibility along the input line. For example, in 
the blunt-body program the enthalpy is determined with pressure and density as 
inputs, whereas in the method of characteristics program the enthalpy is 
determined with pressure and entropy as inputs . In regions where the curve 
fits are poor, particularly near the limits of the tables, the two values for 
enthalpy may be substantially different and may cause difficulties . 

As examples of flow-field calculations for real gases, the pressure dis- 
tribution on the blunt cone-cylinder-flare body studied in the preceding 
example is shown in figure 9 for air and for a mixture of nitrogen and carbon 
dioxide . No unusual effects are noted compared with the earlier perfect gas 
results . 


CONCLUDING REMARKS 


The Ames computer programs for calculating the complete subsonic- 
supersonic flow field around blunt-nosed bodies at zero angle of attack were 
described. The more complex portions of the programs were explained in detail, 
and flow charts for the main programs were presented. The flow charts should 
be helpful in diagnosing difficulties or making minor modifications for 
specific applications. 

A number of example calculations were presented to illustrate the appli- 
cability and accuracy of the programs. It was shown that the use of a quad- 
ratic rather than the usual linear interpolation for entropy improved the 
accuracy of the method of characteristics program. The improvement showed up 
especially in the total pressure near the surface of blunt cones, where large 
entropy gradients develop, and in the integrated mass flow across the shock 
layer for such bodies. Surface-pressure distributions on blunted cones are 
in agreement with published numerical results obtained by somewhat different 
methods. Surface pressures and shock shapes, including the embedded shock, 
for a blunted cone-cylinder-flare, show good agreement with experiment. 
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Finally, to illustrate applicability to flows of real gases in thermodynamic 
equilibrium, surface pressures on a flared body were presented for flight in 
air and in a C 0 2 -N 2 mixture . 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., June 8 , 1965 
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' TABLE II.- LIST OF SUBROUTINES 


BLUNT-BODY PROGRAM 



Subroutine name 

Primary 



and ID no . 

calling 

Subroutine function 


(if any) 

routine 


1 

EXECl(TH070l) 


MAIN 

Reads input cards 

2 

EXEC2(TH0702) 


MAIN 

Reads shock-shape parameters and calcu- 
lates shock shape and slope 


0UT(TH0704) 

1 

r MAIN 

Outputs free-stream conditions and field 

3 

1 

IDERIVT 

data after each step 

b 

SHOCK (TH0705) 


MAIN 

Calculates flow properties behind shock 




r 

wave 

5 

derivs(tho7o6) 

■4 

MAIN 

STEP 

Calculates derivatives in s direction 

6 

FIELDS (TH0707) 


BODYS 

Calculates coefficients for line on 
which output data is desired 

7 

STEP(TH0708) 


MAIN 

Predicts and corrects the flow variables 
for the following step 

8 

B0DYS(TH0709) 


STEP 

Calculates hody location and flow 
proper! ies 

9 

P0LYN(TH0710) 


EXEC2 

Evaluates polynomial and first two deriv- 
atives for given value of argument 

10 

FIELDP(TH071l) 


BODYS 

Interpolates field data to find proper- 
ties on output line 

11 

B0DYS1(IH0712) 


BODYS 

Smooths hody coordinates 

12 

FHIPSI(TH0713) 


STEP 

SHOCK 

Calculates stream function 

13 

terp3(tho7i4) 


FIELDS 

FIELDP 

SONIC 

Interpolates using 3-point Lagrange 
method 

lb 

SONIC (THO 715 ) 


BODYS 

Locates sonic line 

15 

out7(tho7i6) 


MAIN 

Outputs data for hody and along output 
line 

16 

DERIVT(TH0717) 

- 

[main 

(STEP 

Differentiates numerically to obtain 
derivatives in t direction 

17 

FIELDX(TH0718) 


BODYS 

Stores starting data for the method of 
characteristics on tape 

18 

P0LY3(TH0720) 

- 

[MAIN 

jFIELDP 

Evaluates coefficients for second degree 
polynomial passing through three given 
points 

19 

LES2N ( TH072 3 ) 


BODYS 

Evaluates coefficients for least-squares 
fit straight line 

20 

SMOOTH (TH0724) 

4 

fBODYSl 

tDERIVT 

Smooths data hy filtering out high fre- 
quency oscillations 

21 

DERIV1 (TH0725 ) 


fBODYSl 

LDERIVT 

Calculates derivative using 5 points 

22 

DIFF0R(TH0733) 


OUT 

Calculates fourth differences 
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TABLE II.- LIST OF SUBROUTINES 


BLUNT-BODY PROGRAM - Concluded 



Subroutine name 
and ID no. 

(if any) 

Primary 

calling 

routine 

Subroutine function 

23 

ENQ55R(TH0737) 


B0DYS1 

' EXEC1 
SHOCK 

Obtains running integral of equally- 
spaced data 

24 

RGAS 

< 

DERIVS 

BODYS 

FIELDP 

0UT7 

Calculates thermodynamic properties 

25 

SERCH 


RGAS 

Searches for coefficients for real-gas 
properties 

26 

LOCATE 

- 

f FIELDX 
[_RGAS 

Locates tape at specified file 
posit ion 
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TABLE III.- LIST OF SUBROUTINES 


METHOD OF CHARACTERISTICS PROGRAM 



Subroutine name 
and ID no . 

(if any) 

Primary 

calling 

routine 

Subroutine function 

1 

EXEC 


MAIN 

Reads input cards, reads input tapes. 

2 

EXEC2(TH2485) 


MAIN 

initializes variables 
Reads additional input data 

3 

TOP 


MAIN 

Locates nev bow shock point and calcu- 

4 

MID 


MAIN 

lates the new shock angle 

Locates new mesh points and iterates 

5 

B0T(TH2423) 


MAIN 

for solution of equation (10) 
Locates new body point and iterates 

6 

ESHOCK 


MAIN 

for body data 

Keeps track of embedded shocks and 

7 

G SHOCK 


ESHOCK 

adjusts storage locations in P array; 
interpolates for conditions on the 
body upstream of a corner 
Locates general embedded shock point 

8 

C SHOCK 


ESHOCK 

and calculates new shock angle 
Calculates shock angle at a corner, and 

9 

EXPAN 


ESHOCK 

at the first mesh point away from the 
corner 

Calculates additional points for an 

10 

SHOCK 


'top 

GSHOCK 

expansion corner, and adjusts storage 
locations in P array 
Calculates the shock jump conditions 
given the shock angle and upstream 

11 

PM2 


.C SHOCK 
EXPAN 

conditions 

Calculates Prandtl-Meyer flow given the 

12 

13 

RGAS 

ROOTB ( TH2 42 6 ) 


BOT 

expansion angle and upstream conditions 
Calculates gas properties; reads RGAS 
tape 

Locates intersection of right-running 

l4 

C0N(TH24ll) 


MID 

characteristic with the body 
Calculates averages of the coefficients 

15 

DATA ( TH2 42 5 ) 


EXEC 

of equation (10) 

Reads starting flow-field data if 

16 

17 

TPRES(TH2429) 

ISENC 

-< 

f TPRES 
1PM2 

specified on cards 
Calculates total pressure 
Calculates one -dimensional isentropic 
flow between given velocities 

l8 

ESPACl(TH2409) 


EXEC 

Prepares data for equal spacing 

19 

ESPACE (TH2412 ) 


ESPAC1 

Equally spaces data with respect to a 

20 

21 

NTERP 

TERP3(TH2405) 



given variable 
Interpolating routine 
Interpolating routine 
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TABLE III.- LIST OF SUBROUTINES 


METHOD OF CHARACTERISTICS PROGRAM - Concluded 



Subroutine name 
and ID no . 

(if any) 

Primary 

calling 

routine 

Subroutine function 

22 

23 

24 

25 

TERP4(TH24o6) 

POLY3(TH24o8) 

HERM 

ENQ55R(TH0737) 

1 


Interpolating routines 
Integration routine 

26 

shockp(th2487) 

MAIN 

Interpolates for data at prescribed 
points along the bow shock 

27 

B0DYP(TH2488) 

MAIN 

Interpolates for data at prescribed 
points along the body 

28 

SEARCH (TH2489) 

MAIN 

Interpolates for data along a charac- 
teristic line 

29 

SERCHl(TH24l6) 

BODIP 

Calculates equations of body normal 
probe lines 

30 

PRINTF(TH2486) 

MAIN 

Stores data along body or axis normals; 
stores last probe on tape 

31 

SERCH 

RGAS 

Scans RGAS table for pertinent data 

32 

PLOTS 

MAIN 

Dummy subroutine - can be used to 
write a plot tape using data in P 
array 

33 

LOCATE 

(exec 

\PRINTE 

Locates specified file positions on 
data storage tape 
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Figure 2.- Flow chart for characteristics program 
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(c) Moq = 10 
Figure 5*- Concluded. 
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O Experiment, ref.- 18 
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